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We study the relativistic heat equation in one space dimension. We prove a local regularity 
result when the initial datum is locally Lipschitz in its support. We propose a numerical scheme 
that captures the known features of the solutions and allows for analysing further properties of their 
qualitative behavior. 
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53 ! 1 Introduction 

In this work, we explore both analytically and numerically the implications of a new strategy to study 
flux-dominated nonlinear diffusions in one dimension. To be more precise, we consider the so-called 
1 relativistic heat equation (RHE) 

en" 
m 




Ut = u \ - UUx I , x G R, t > 0. (1.1) 



■ introduced by Rosenau in [37] and, later on, by Brenier in p3| based on optimal transportation ideas. 
The name of RHE comes from the fact that (II. ip converges as c — > oo to the heat equation both formally 
and rigorously [20], while the flux in ([Lip , understood as a conservation law, is bounded by the speed 
of light c whenever the solution is positive. 

Many other models of nonlinear degenerate parabolic equations with flux saturation as the gradient 
| becomes unbounded have been proposed by Rosenau and his coworkers |19[ 137]. and Bertsch and Dal 
Passo [121 [26] . Notice also [36] for the presence of flux limited diffusion equations in the context of 
radiation hydrodynamics. 

The general class of flux limited diffusion equations and the properties of the relativistic heat 
equation have been studied in a series of papers [5j HI [6j [21] . An existence and uniqueness theory of 
entropy solutions for the Cauchy problem associated to the quasi-linear parabolic equation 

Flit 

— = &wh(u,Du), (1.2) 

was developed in [5H3]. Here, the flux function is given by b(z,£) = V^/(z,^) and f : Mx M N — > JR + 
is a convex function with linear growth as ||£|| — > oo, such that V^/(^, £) G C(M x M N ) satisfying 
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other additional technical assumptions. In particular, the relativistic heat equation (jl.ip satisfies these 
assumptions, and other models considered in [37]. To avoid the difficulty of the lack of a-priori estimates 
that ensure the compactness in time of solutions of (|1.2|) . the existence problem was approached using 
Crandall-Liggett's theorem [24]. For that, we first considered the associated elliptic problem and we 
defined a notion of entropy solution for which we developed a well-posedness theory. The notion of 
entropy solution permits to prove a uniqueness result using Kruzhkov's doubling variables technique 
|30j [T5] . This technique was suitably adapted to work with functions whose truncatures are of bounded 
variation [U 0], which is the natural functional setting for (II. 2D and its associated elliptic equation. 

The evolution of the support of solutions of the relativistic heat equation (jl.ip was studied in [6]. 
By constructing sub- and super-solutions which are fronts evolving at speed c and using a comparison 
principle between entropy solutions and sub- and super-solutions, it was proved in [6] that the support 
of solutions evolves at speed c. Moreover, the existence of solutions which have discontinuity fronts 
moving at the speed c was again shown using the comparison principle with sub-solutions. This implies, 
in particular, that the maximal regularity in time that one can expect for general solutions of (jl.ip is 
that u G BV([t, T] x M n ) for any < r < T. That this happens for a general class of initial conditions 
was proved in [8] and later extended in [21] . This lack of regularity is at the origin of the notion of 
entropy solutions for this type of equations. The only regularity result for smooth initial conditions 
was proved in [20] and it guarantees that V In u is bounded whenever initially is. But the study of the 
local regularity of solutions of (II. ip is still an open question. One of the purposes of this paper is to 
address this problem for the Cauchy problem associated to (jl.ip in one space dimension with compactly 
supported bounded probability densities as initial data. 

Assuming that the initial data in non-negative, we can easily change variables to observe that u(t, x) 
is a solution of (jl.ip if and only if u(t, x) = u(-^t, ^x) is a solution of 



Thus, without loss of generality we may assume that v = c = 1, and for simplicity we shall assume 
it in the rest unless explicitly stated. Notice also that if u(t) is a solution corresponding to uq, then 
Xu(t) is a solution corresponding to Xuq, A > 0. Thus, without loss of generality we assume that 
H^o 111 = || u(t) ||i = 1 for any t > 0, and reduce our evolution to probability densities. In this paragraph, 
the term solution refers to entropy solution for which the well-posedness theory was developed and for 
which a summary of its concept is reminded to the reader in the Appendix. 

The local regularity of entropy solutions to (jl.3p will be done by a change variables, writing (|1.3|) in 
terms of its inverse distribution function. This change of variables has its origin in using mass transport 
techniques to study diffusion equations |18j [T3] . It is known p3] that equation (jl.ip has the structure of 
a gradient flow of a certain functional (the physical entropy) with respect to some transport distance. 
This structure was already used to give well-posedness results to (jl.ip in [35] • Nonlinear diffusions have 
received lots of attention from optimal transport theory viewpoint starting from the seminal works 



Transport distances between probability measures in one dimension are much easier to compute since 
they can be written in terms of distribution functions and their generalized inverses (pseudo- inverse), 
the so-called Hoeffding-Frechet Lemma |39j Section 2.2]. This result led to the following change of 
variables based on the distribution function F associated to the probability measure u, defined as 




(1.3) 



[291 [33]- 
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We formally consider its inverse (p denned on the mass variable r\ G (0, 1) that verifies 



F(t,tp(t,V))=V, ??e(0,i)- 

After straightforward computations assuming that all involved functions are well-defined and smooth, 
one obtains the equation 

= 'jm (1.4) 

for the inverse distribution function (p. This change of variables has first been used for nonlinear diffu- 
sions in [18] to show contractivity properties of transport distances for porous-medium like equations. 
It is worthy to remark that an implicit Euler discretization of (11. 4j) is equivalent to the variational JKO 
scheme whose convergence is proved in [35] for (jl.ip under certain assumptions. Numerical schemes 
to solve the equation for the pseudo-inverse function in the case of the porous medium equation were 
analysed in |28j . This Lagrangian approach was generalized to several dimensions in [16] in order to 
propose numerical schemes for equations with gradient flow structure in optimal transport theory and 
general quasilinear problems in divergence form. 

In Section [21 we will first take advantage of this change of variables to prove the following regularity 
result: 

Theorem 1.1. Let u$ G L°°(]R) with uq(x) > k > for x G [a, b], and uo(x) = for x £" [a, b\. Assume 
that uq G VF 1,00 ([a, b]). Let u(t,x) be the entropy solution of (jl.ip with u(0) = uo, ||«o||i = 1 • Then 
u G C{[0,T},L 1 {M N )) satisfies: 

(i) u(t,x) > K,(t) > for any x G (a — ct,b + ct) and any t > 0, u(t, x) =0, i ^ [a — ct, b + ct\, 
t€(0,T), 

(ii) u(t) G BV(¥L), u(t) G W 1 ' 1 {a - ct,b + ct) for almost any t G (0,T), and u(t) is smooth inside its 
support, 

(Hi) if uq G W 2,1 (a,b), then Ut is a Radon measure in (0, T) x M. 

We emphasize that the new parts of this result with respect to the literature discussed above refer 
to the regularity stated on points (ii) and (iii). This result implies that sharp corners on the support of 
the initial data are immediately smoothed out by the evolution of the RHE. This result will be extended 
in Section [3j in particular, we cover the case where the initial condition uo vanishes at the boundary of 
its support. 

In Section^ we will propose an adaption of the numerical scheme in [16] based on equation (|1.4p 
with suitable boundary conditions that fully captures the demonstrated behavior of the solutions of 
the RHE. Moreover, we will show different numerical tests in situations where the theory has not 
been developed yet. For instance, we numerically study the conditions for the formation or not of 
discontinuities on the bulk of the solutions for RHE and its porous medium counterparts 



u t 



y/ U 2 + { Ux y 



with m > 1 and their long-time asymptotic behaviour. Finally, we include in Appendix A some basic 
material to describe the notion of entropy solutions for (jl.3p for the sake of completeness. 
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2 Regularity of Solutions 

As proved in [5], there exists a unique entropy solution of the Cauchy problem for (jl.3j) for any uq G 
L 1 (iR)nL°°(JR), no > 0, see Appendix for the full notion of solution. Moreover if uq has compact support 
in M and is locally bounded away from zero in any interior point of its support, then supp(u(i)) = 
supp(iio) © B(0,t) [6]. The rest of this Section is devoted to the proof of the regularity statements (ii) 
and (hi). 

Let us recall that the entropy condition on the jump set of u can be expressed by saying that the 
profile of u is vertical at those points. Since the support of u(t) is (a — t,b + 1), and u(t) > n(t) > in 
(a — t, b + 1) for any t > [6], there is a jump at the points x = a — t,b + t and we have [21] 



= (t,o-t) = l, . n U * (t,b + t) = -l. (2.1) 



Let us consider the change of variables discussed in the introduction and define the function ip(t, rj) by 
the relation 

r<p(t,rj) 

u(t,x)dx = r], r) G (0,1). (2.2) 

la-t 

Proceeding formally, assuming that the function is smooth inside its support and differentiating with 
respect to r\ we obtain 

u(t,x)ip v = 1, for x = tp(t,r)). 
Differentiating with respect to t we have 

u(t,x)(f t + u(t,a — t) + / Ut(t, x)dx = 0. 

J a—t 

Taking into account the boundary conditions (|2.ip |21] , one has 

<p(t,r)) r<p(t,v) 



a—t 



u t (t,r)dr = I I / da;= J^ =(t,o;)-M(t,a-t), 



hence 

n(t, x)v9 t = x ' (t, x) for x = <p(t, rj) . 

^u 2 + (u x ) 2 

Then the equation satisfied by <p is 



2.1 Regularity result in mass variables 

Now, let us consider the change of variables v = ip v . The equation satisfied by v is 

v t = 1—=^==) t>0, 16(0,1). (2.3) 

where we have written x instead of rj. This will done through this subsection for convenience. 
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The initial condition vq is determined from the initial condition uq. We assume that uq G L°°(iR), 
uq > k, and wo G VF 1,oc, ([a, b]). Since the relation between uq and v is determined by vo(n) = , 

then ai := tt-V- < i>o < ^ := a 2 . We also have u G VF 1,oo (0, 1). Note that 

1 /-b 
vo(r]) dr] = / dx = b — a. 



o Ja 

If we denote by z/ the outer unit normal to (0, 1), that is z/(0) = —1 and z/(l) = 1, the natural boundary 
conditions for (|2.3[) are 

i = m v = l at xe 5(0,1), (2.4) 

with 9(0, 1) = {0, 1}. The first step toward Theorem II .11 is to show a regularity result for the Cauchy 
problem ([Ol - dS^l) . 

Theorem 2.1. Assume that vo G M^ 1,oo (0, 1), vo > a\ > 0. Then there exists a smooth solution of 
(|2.3p m (0,T) x (0,1) u>it/\ u(0,x) = vq{x) and satisfying the boundary conditions (|2.4p fin a weak 
sense). 

Proof. To prove this claim, we consider the following approximated Cauchy problem 

v t = l-—^==\ +ev xx t€ (0,r), xe (0,1), (2.5) 

( . , 2 + eu x U = l-e 1 / 3 , t€(0,T),i€S(0,l), (2.6) 

where e > 0. The proof is divided in several Steps. In Steps 1 to 3 we prove some formal estimates 
that are also useful to state the existence of solutions of (|2.5p - (|2.6p in Step 4. For simplicity we write 



&{z, = - z > 0, £ G M. 

Let us observe that 

a(z,£)£>|£|-z 2 . (2.7) 

Step 1. LP bounds on v for p G [l,oo). Let us first consider the evolution of the L 1 norm. For that we 
integrate (|2.5p on (0, 1). We have 

d 

and thus, 



l 

^ / v(t,x)dx = (a(v,v x ) + ev x )(l)-(a(v,v x ) + ev x )(0)=2(l-e 1 / 3 ), 
dt J 

[ v(t,x)dx= [ v (x)dx + 2(l - e 1/3 )t. (2i 
Jo Jo 



Given 1 < p < oo, we have 
1 d 



f 1 v p+1 (t,x)dx+ f 1 a(v,v x )(v p ) x dx + ep C v p ~ l {v x f dx = (1 - e 1 / 3 ) [ 

Jo Jo Jo Jd( 



p + Ldt Jo Jo Jo Jd(0,l) 

< f v p dx+ [ \(v p ) x \dx, 
Jo Jo 
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where the inequality 

v p (0) + v p (l)= [ v p < [ v p dx+ [ \{v p ) x \dx 

Jd{0,l) JO JO 

holds in one dimension. Using (j2.7[) we have 

a(v,v x )(v p ) x > [ \(v p ) x \-p f v p+1 , 
Jo Jo 



hence 

1 d 



v p+1 (t,x)dx + ep v p ~ l {v x ) 2 dx< \ v p dx + p v p+1 dx 



p+ldt 

Using this recurrence relation, by Gronwall's inequality, we obtain that 



f v{t,x) p dx <C(T,p) Vt € (0,T),Vp € [l,oo) 
Jo 



and that 



vP^ivxfdxdt < C(T,p), Vp€[l,oo), (2.9) 
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where the constant C(T,p) does not depend on e. 



Step 2. L°° bounds above and below on v independent of e. Let us construct a supersolution to the 
Cauchy problem (|23j) - pi)]) . Let V{t,x) = B(t) - V / e^ 

+ x(l — x) with B smooth and increasing. 

Take B(0) such that 

V(0,x) = B(0) - y^ 2 / 3 +x(l - x) > v (x). 



We compute 



V t = B'(t), 

_ (x - 1/2) _ e 2 /3 + 1/4 

"x — ; — ^-r-- 1 "a 



y / € 2/3 + x{1 _ x y (£2/3+x(l-x)) 3 /2' 
Ur (S - 1/2) 



a^K) 



(U 4 + ^2)1/2 " D (t } 
M 2 



x 



where D(t, x) = (V(t, x) 4 (e 2 / 3 + x(l - x)) + (x - 1/2) 2 ) . Note that D(t, x) is a smooth and strictly 
positive function in [0, 1]. Moreover, since B is increasing, D > (U(0,x) 4 (e3 + x(l — re)) + (x — ^) 2 )^. 
Thus |a(V, V^)a;| < C for a constant C that can be taken independent of e and t G [0, T]. Thus, a direct 
computation shows that 

,2/3 _|_ 1 /A 1 ^ 

a(U, V x ) x + eU^ < C + e < C + e 2 / 3 + -<C + - = C, 

where C does not depend on e £ (0, 1]. Take B'{t) > C, for instance B(t) = B(0) + Ct. Let us prove 
that given T > 0, for e > small enough V(t, x) satisfies 

(a(U,y) + eU>> 1-e 1 / 3 , 

for t £ [0,T], hence V(t, x) is a supersolution of the Cauchy problem (|2.5p - (j2.6p in [0, T]. Indeed, since 

D(t,0) = f (fi(t) - e 1 ^ 2 / 3 + 1/4 N 1 
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we have at x = 

(a(V, 14) + e y,),U=o = ^ + = — —n- 2 + ^ 2/3 > 1 - . 1/3 

D(t,0) e V3 (l + 4( J B(t)-e 1 /3)4 e 2/3) 1 /2 2 

for e > small enough, and analogously at x = 1. Since V(t, x) is a supersolution for the Cauchy 
problem (|2.5p - (|2.6p . by the classical comparison principle we get v < V in [0,T] x [0, 1], and thus there 
exists M > depending only on uq and T such that v(t, x) < M in (i, x) G [0, T] x [0, 1]. 

Let us finally observe that v > a±. Indeed, v = a\ is a subsolution for the Cauchy problem 
(|2.5p ~ (|2.6p and by the comparison principle in its weak version, we deduce that v > a±. 



Step 3. L p bounds on v x independent of e. Putting together the estimates in Step 2 and (|2.9p . we 
deduce that 

\(v p ) x \dxdt<C(T,p), 

Jo Jo 
for any p E [1, oo). 



Step 4- Existence of smooth solutions for the Cauchy problem (I2.5p - (l2.6p . The existence of solutions of 
(|2.5p - (|2.6p follows from classical results in [31] and |32| Theorem 13.24]. We note that thanks to the a 
priori bounds stated above we could use the flux 

a M (v,v x ) 



^/inf (H,M) 4 + n 



X 



so that the assumptions of the existence theorems in |31j and \32\ Theorem 13.24] hold. Finally, observe 
that we need to assume a compatibility condition on vo so that vq satisfies (|2.6p . If does not satisfy 
(|2.6p . we modify it to define a function Vo, e G VF 1,oo (0, 1) satisfying (|2.6p . This modification is only 
done in a neighborhood of x G <9(0, 1) which vanishes as e — > 0+, so that VQ )t is locally Lipschitz inside 
(0, 1) with bounds independent of e. Finally, we observe that this modification can be done in such a 
way that 

sup eH^oezlloo < oo. ( 2 -10) 

£6(0,1] 

Although we omit the details of the construction, let us check that (|2.10p is compatible with (|2.6|) . 
For that, notice that we can take vo €X = A(e)e~ a with a = | and A(e) = -^vq 6 (0) 2 + C^e 1 / 3 ). Indeed 
substituting this expression in (12.6p . we have 

Me)/t a A(e) _ 1 ^ 



V«o £ (0) 4 + A(e) 2 /e 2a 



An asymptotic expansion shows A(e) = -^u 0e (0) 2 + Oie 1 / 3 ), and thus (pTOl) is compatible with (pToD . 

Let v e be the solution of the Cauchy problem (l2.5p -( i2T6j) . Then v t has first derivatives Holder 
continuous up to the boundary and for g = v exx , v e t, we have 



sup 

x^y 



min(d((x, t),V),d((y, a), V)) l ~ 8 \^9^)\ 

[\x — y\ z + \s — t\) a l' L 



for some a, 5 > 0, where V is the parabolic boundary of (0, 1) x (0, T), that is [0, 1] x {0}U{0, 1} x (0, T), 
and d(-,V) denotes the distance to V . On the other hand, by the interior regularity theorem |3H 
Chapter V, Theorem 3.1], the solution is infinitely smooth in the interior of the domain. At this point 
the smoothness bounds depend on e. 
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Step 5. A local Lipschitz bound on v e uniform on e. For simplicity of notation, let us write v instead 
of v e . Let w = |f x | 2 (/> 2 where <\> > is smooth with compact support [0,6] C (0,1). This Step is a 
consequence of the following inequality 

w t < A(t, x)w xx + B(t, x)w x + Cw + f(t, x), (2.11) 

where A, B are smooth functions, C = (12 + |), and < / = P(v,ip,<p x )\<p x \ + ^e(f) 2 v 2 , where P is a 
polynomial in v of degree 3. Assume for the moment that the last term ellu^^Hoo £ L°°(0,T). Using 
Step 2, this implies that / G L°°([0, T] x [0,1]). Thus we may replace / by ||/(i)||oo- The change of 
variables 

w(t, x) = e~ ct w(t, x) — / f(s)ds 



permits to write (|2.1ip as wt < A(t, x)w xx +B(t, x)w x . Then, using the maximum principle, this implies 

sup ||u)(£)||oo < ^(O)^ 

te[o,T] 

hence we get 

sup |K*)||oo < C(T,0, IKOJUoo). 
te[o,T] 

Let us now prove the claim (|2.1ip . We first compute 



( z 4 + £2)3/2' "V >sy ^4 + £2)3/2 (^4 + £2)5/2 ' 

a g(^£) = 7-4- ^ 2 \3/2 ' a €z(^£) = / 4 , ^5/2 1 and a «(^'0 



( z 4 + £2)3/2' t*V ( z 4 + £2)5/2 ' «V (^4 + £2)5/2 ' 

We also compute w x = 2(/></> x ?; 2 + 2^ 2 t! ;r t> : r:r and w xx = (24> 2 X +2<p<p xx )v 2 + 8(p(p x v x v xx +2(p 2 v 2 x +2<j) 2 v x v xxx . 
Differentiating (j2.5[) with respect to x and multiplying by <ft 2 we obtain 

-w t = a zz vl(f) 2 + 2a^ z v 2 x v xx ^ 2 + a^v x v 2 xx (j) 2 + a z v x v xx (f> 2 + a^v x v xxx (p 2 + ev x v xxx ct> 2 . 
Now, we get 

3 2 _ 6^ 2 t;^ 2 \2v % v x 4> 2 

a ** V *<P ~ - {v 4 + v 2J3/2 + {v 4 +v 2 r /2 $ ^ , 

2a^ z vlv xx (j) 2 = a^ z v x w x - 2a£ z vl<jxj) x < a^ z v x w x + 12v 3 (p\4> x \ , 
a^v x vl x 4> 2 = ^a^v xx w x - a^v xx v 2 x ^ x = ^a^v xx w x - X, 



and 



where X = a^v xx v x (fxj) x . Similarly, we obtain 



and 



a- z v x v xx (j) 2 = ^a z w x - a z vl<p(j) x < ^a z w x + 2-y 3 ^, s 



^v x v xxx (f) 2 = ^a^w xx - a ? (</> 2 . + 4>4> xx )vl - Aa^v xx v x 4>4> x - a^v 2 ^ 2 
< ^w xx + v 2 ((f>l + </)\4> xx \) - Y - a^cf) 2 , 
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where Y = 4a^v xx v x (p4> x - Direct estimates show that 

\Y\ < ^u,i,,:r + < \^v 2 xx 4> 2 + 8v 2 <g 

and 2 

I VI z 1 2 /2 , a I? 4 j.2 ^ ^ 2 j.2 , ^ 2/2 
1^1 < ^H V xx<^ + ^ V x^x < ^V xx <j) + -V <j> x . 

Finally, let us compute the term 

V x Vxxx4> 2 = ^jWxx ~ (4>l + 4>4>x)vl ~ 4:(f)(f> x V x V xx - (j) 2 V 2 xx 



< ^w xx - 4> 2 x v 2 x + -4> 2 v 2 x + -4> 2 x v 2 x + 4<j) 2 x v 2 x + 4> 2 v 2 xx - 4> 2 v 2 xa 



ill 1/22 1 

-w xx - <p x v x + -<j> v x + -< 

1 1 7 j2 2 

= ~w xx + 2 W + 29x v x- 
Putting all together, we get the desired claim (|2.1ip 

^w t < ^(a^ + ejwxx + (a.£ Z v x + ^a.^v xx + ^a^j w x + (l2 + w + P(v, <j>, <j> x )\<i>x\ + ^4> 2 x v 2 x , (2.12) 

where P is a polynomial of degree 3 in v. 
Now, we have to show that ell 

(t) || oo £ (0, r) . Let us first exploit the boundary condition in 
(|2.6p . Multiplying it by v x and using (|2.7p . we get 

\v x \ - v 2 < a(v, Uj,.)^ = " . + ev 2 x = (1 - e 1 / 3 )^ , 

and thus we get that ev 2 < v 2 on 9(0, 1). Moreover, using Step 2 we finally deduce that 

ev 2 x (t) < sup (K(i, 0)|, 1)|) < M, on 9(0, 1) . (2.13) 

te[o,T] 

Taking = 1 in (|2.12p , we obtain 



^w t < ^(a ? + e)w xx + (a (z v x + ^a^u^ + ia^ u^. + (12 + |; 



that together with (|2. 13|) and the maximum principle, imply that 

e|Mt) 2 ||oc < C, (2.14) 

for some constant C that depends on the bound (|2.10p . and thus independent of e. 

Summarizing, now the term |e||0 2 ^(t)||oo £ L°°(0,T) with bounds independent of e. Again, Step 
2 implies that ||/(i)||oo < \\P(v(t), 4>, 4>x)\<Px\\\oo + ll|e<#^(*)||oo € L°°(0,T) with bounds independent 
of e. Then the argument given above shows that there are local Lipschitz bounds on v e uniform in e. 

Step 6. Interior regularity of higher order derivatives uniform in e. Thanks to the smoothness results 
stated in Step 4 and the local uniform bounds on the gradient in Step 5, the classical interior regularity 
results in |31} Chapter V, Theorem 3.1] shows uniform (in e) interior bounds for any space and time 
derivative of v f . 
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Step 7. Passing to the limit as e — > + . Letting e — > + is not completely obvious due to the boundary 
condition (I2.4p . Another difficulty stems from the fact that we do not know if v e t are Radon measures 
with uniform bounds in e. This means that the notion of normal boundary trace has to be considered 
in a weak sense as considered in [2] (see also (3J Section 5.6] or [9]). Thus, we only sketch the proof of 
this result. Let us first prove that the interior regularity bounds on v e permit to pass to the limit and 
obtain a solution v of 



Let 



tH= ^ ■ inP'((0,r)x(0,l)). 



£ e := v et = [ r== = + ev €X ) and a e = j== = + ev 



Estimate (|2.14p implies that a e are uniformly bounded independently of e. Then by extracting a 
subsequence, we may assume that a e — a G L°°((0,T) x (0,1)) weakly*. On the other hand, the 
interior regularity bounds on v e ensure that a = , =5 =. By passing to the limit as e — > 0, we have 

V « +(«!!! ) 2 

v t = a x in ©'((0, T) x (0, 1)). Finally, if we take ip G C^QO, T] x [0, 1]) with y>(0) = tp(T) = 0, multiply 
(|2.-5[) by ip and integrate by parts, we obtain 



[ [ v e ip t dxdt= [ [ & e ip x alxdt - 2(1 - e 1/3 )T. 
Jo Jo Jo Jo 



Letting e — > + , we obtain 



T fl pT i-l 

vipt dxdt = I / sap x dxdt — 2T. 
Jo Jo 



This is a weak form of the boundary condition (|2,4p . The correct notion of weak trace is much more 
technical and is described in [3]. Using Lemma 5.7 in [9] one can directly obtain that v satisfies (j2.4|) 
in this generalized sense. Since we do not need this result here, we skip the details that would need 
several technical definitions to be fully explained. □ 

Remark 2.2. Note that we can apply Step 5 to the smooth solution obtained in Theorem 12.11 to the 
Cauchy problem (|2.3p - (|2.4p . In this case ||/||oo < ll-P( v > <^> ) I I ||oo and we obtain a local Lipschitz 
bound for v(t,x) which only depends on local uniform bounds of v(t,x) and on the local Lipschitz 
bound of vq(x). 

Remark 2.3. In Section 12.21 we will give sufficient conditions on uq that imply that vt is a Radon 
measure. In that case, the notion of weak trace a • v can be found in |21l 123] . 

Remark 2.4. We could define the notion on entropy solutions of equation (|2.3p with boundary con- 
dition (|2.4p and prove that the solution constructed is indeed an entropy solution of it. We will not 
pursue this here. 

2.2 Getting an entropy solution of (11. 3p from ( 12. 3 j) 

Here, we use several notations and definitions that are introduced in the Appendix to which we refer 
for details. In this Section, we come back to the notation v(t,rj) instead of v(t,x), r\ G (0, 1). Recall 
that by passing to the limit as e — > + we have found a solution v of 

v t = ( / ^ 7 ~ ] inP'((0,r)x(0,l)), (2.15) 
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for any T > 0. Thus, let v(t, rj) be the solution of (|2.15p constructed in Theorem 12.11 which satisfies 
[a(i, rj) • v\ = l for rj = 0, 1 and a.e. for £ £ (0, T) in a weak sense. As we shall see, we do not need this 
here, we only need a weaker form of the boundary condition as expressed in (|2. 17j) below. 

In the next Lemma we construct an entropy solution of (jl.ip from a solution v(t,rj) of (I2.15p . To 
prepare its statement, let no £ L°°(R) with Uq{x) > k > for x £ [a, 6], and no(x) = for x £" [a, b]. 
Assume that no £ W 1,OQ ([a, b]). Let vq(t]) = ^4^] > ^ ^ (0j 1)j where x = (p(0,r)) is such that 

r<p(o,v) 

/ no(x) dx = rj . 

J a 

Let n(t, x) be defined in [a — f , b + t] by 

1 r _ _ 

u(t,x) = — -, where x = tp(t, rf) = a — t + / v(t,r])dr). (2-16) 

v{t,V) Jo 

By (J23D, we have 

-l 

v(t,rj)drj = b-a + 2t, (2.17) 



and x = </?(i, 77) £ [a — t, 6 + t] when r? varies in [0, 1]. Note that 



u(t, x)dx = rj, rj £ (0, 1). 



a—t 



We define u(t,x) = 0, x £" [a - 1, b + t], t £ (0,T). Notice that u(t,x) > > for any x £ (a — t, b+t) 
and any t > 0. 

The statement (ii) in Theorem 11.11 follows from next Proposition. 

Proposition 2.5. Given u defined by (|2.16p where v is a solution given by Theorem 12. 1 L T/ien u £ 
C([0, T], L 1 (JR)), u(0) = no, and satisfies 

(i) u{t) £ BV(M), u(t) £ W 1 ' 1 ^ - t,b + t) for almost any t £ (0, T), and u(t) is smooth inside its 
support, 

(ii) u t = z x m V'((0,T) x R), where z(t) = "S" z(f) m2 , 

(Hi) u(t,x) is the entropy solution of (jl.3p with initial data uq in (0,T). 

Proof, (i) Since v is bounded and bounded away from zero from Step 2 in Theorem I2.1| then u is 
bounded and bounded away from zero in its support. The smoothness properties of v prove that 
n £ C([0, T], L l (R)), n(0) = no, and u(t) is smooth inside its support. By Step 3 from Theorem 12. 1| 
we have that u(t) £ W 1 ' 1 (a -t,b + t) for almost any t £ (0,T). This implies that n(i) £ BV(R) for 
almost any i £ (0, T). From the change of variables f)2. 16|) we have that 

(2.18) 



(m) For simplicity, let us write Qt = (0, T) x J?, and f2(t) = (a — t, b + i). Since 

Du(t) = u x X m - u'ftuWl-dSlit), 
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we have that u G Lj oc W (0,T; BV(M)). We have denoted by u l (t) the trace of u\ n ^ on dU(t). Note 
that it coincides with u + (t). Let us prove that 

u t = z x in V((0,T) x M). (2.19) 

Let <j) G V(Q T ). Let 77) = <p(t, rj)), rj G [0, 1]. Then t = <fo(i, rj)) + x (t, <p(t, 77))^ and 

/ u<ptdxdt = — I I u<j>t dxdt = — J j -(4> t — (f> x (t,<p(t,rj))(p t )v drjdt 
o Jm Jo Jn(t) Jo Jo v 

T pi pT pi 

/ (<p t - (f) x (t,(f(t, rj))(p t ) drjdt = / / (f> x (t,<p(t,ri))<p t drjdt 
o Jo Jo Jo 

T f 1 v f T f uu x 

l <p x (t,tp(t,rj)) 2 drjdt = - / x 6 x (t,x) dxdt 

Jo J v i + v 2 Jo Jn(t) V n + u x 

T f 

/ z(p x dxdt , 

o J Ml 

where (pTL8|) was used. Thus (f2T9]) holds. 

(Hi) To prove that u is an entropy solution of (|1.3p . we have to prove that 

hs(u, DT(u))(f> dxdt + / hr(u, DS(u))4> dxdt 
Qt JQt 

< / J T s(u)(p t dxdt- / z(t,x) -V(f>(t)T(u(t))S(u(t))dxdt, (2.20) 

holds for any any T,S G 7~ + and any G 2?((0, T) x JR), <p(t,x) = rj(t)p(x). As in [61 Proposition 1], 
we have 

(h s (u(t),DT(u(t)))) s = \D j J S KT>(u(t))\ = JsRT'WWU^d^t) (2.21) 

and 

(h T (u(t),DS(u(t)))y = \Dij TRS ,(u(t))\ = jTRS'(u l (t))n ^dn(t), (2.22) 

where R(r) =r,r G M. Thus, by (f!T2Tj) and (f!T2"2"|) . we get 

(^(tJ.DTKt))))' + (fc T (u(*) > Z>S(u(t))))* = (J SJH v («<(*)) + J Tft s' («*(*))) U° \_dSl(t) 

= (TSRtfit)) - Jrstu^t))) H°^dn(t). (2.23) 

On the other hand, it is easy to prove that 

(h s (u,DT(u))) ac 0dxdt + [ (h T (u,DS(u))) ac ^dxdt 

Jq t 

[ z(t,x) ■ [T(u(t,x))S(u(t,x))] x 4>(t)dxdt. (2.24) 

'0 Jfl(t) 

Adding (pT23|) and (pT24"|) . we obtain 

<j)h s (u(t),DT(u(t)))dxdt+ / (j)h T (u(t),DS(u(t)))dxdt 

Jq t 

[ (TSR(u\t)) - Jrs&tt))) 4>(t) d%° dt 

JdQ{t) 

+ [ [ z(t,x) ■ [T(u(t,x))S(u(t,x))] x (j)(t)dxdt. (2.25) 
Jo Ju(t) 



12 



To simplify the subsequent notation let us denote p{u) = T(u)S(u) = J'(u) and J(u) = Jrs(it)- 
Let us now prove that 

f (p(u\t))u\t)- J(u i {t)))<j){t)dU dt+ I [ z-[p(u)] x 4>dxdt 
o Jdn(t) Jo Ja(t) 

< / J(u)<pt dxdt — / / <p x zp(u) dxdt. (2.26) 
Jq t Jo Ju(t) 

The main technical difficulty comes from the fact that we do not know that ut = z x is a Radon 
measure. We circumvent this difficulty by using instead discrete derivatives. Let us denote 

1 1 

A+w(t) = -{w{t + r) - w{t)), A~u;(i) = -(w(t) - w(t - r)). 

r r 

Then, we can obtain 

/ up(u)4> A~Xq(£) dxdt = / / (up(u)4>) dxdt 

o Jm Jo Jn(t) 

f A+u(t)p(u(t + T))</>(t + T)dxdt+ f f u(t)A+(p(u)(f))(t)dxdt 
o Jn(t) Jo Jn(t) 

> f [ A+J(u)(t)(j)(t + T)dxdt+ [ [ u{t)A+{p{u)4>)(t)dxdt 
Jo Jn(t) Jo Jn(t) 

J{u){t)A~[<p(t + T)]dxdt- [ f J(u)(t)(f>(t)A^X n(t) dxdt 
'o Jn(t) Jo Jm 

+ f [ u(t)A+(p(u)<f>)(t)dxdt 
Jo Jn(t) 

which is a discrete version of (|2.26p . Note that we have used the inequality A^u(t)p(u(t + r)) > A+J(u) 
which is a consequence of the convexity of J. By letting r — > + , we need to show that 

(u(t)p(u(t)) - J(u(t)))<t>(t)A~X m dxdt -> [ T [ (p(^)K(i) - J(u l (t))) dn° dt, (2.27) 
q t Jo JdQ(t) 

[ J(u)(t)A~[(p(t + T)}dxdt-> [ [ J(u)(t)<f)t(t) dxdt, (2.28) 
o Jn(t) Jo Jn(t) 

and 

fT r rT 



u(t)A+(p{u)4))(t)dxdt / / (p(u)(j)) x zdxdt. (2.29) 
o ^o(t) Vn(t) 

This will result in $2?26\\ . The limit (pT2T)) follows since u(t) G BV(1R) a.e. in t, u G L^(0, T; BV(1R)) 
(hence ||t/a;(i)|| G L 1 (0, T)) and the trace functions u(t,a — t), u(t,b + t) are integrable in [0, T]. The 
second limit (|2.28p follows easily. To prove (|2.29p . for any r > let 

1 [ t+T 

tj) T (t, x) := — / <p(s, x)p(u(s, x)) ds, 

T Jt 

and observe that 

At(p(u)4>)(t,x) = -^r(t,x). 
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Observe also that 

-[r (t, <p(t, 7]))} = -^ r (t, <p(t, ri)) + i>i(t, V {t, v ))<p t (t, v ). 

Then, as r — > + 

T r f T f d 

u(t)A+(p(u)(f))(t)dxdt = I I u(t)—ip T (t,x)dxdt 



JQ(t) J JCl(t) 



Jo 



T r l v 

^l(t,(p{t,r])) 1 drjdt 
o Jo Jv A + v% 

T rl rT 

UU, 



ipl(t,x)—z dxdt -» / / (p(u)cf>) x zdxdt. 

JO yJU~ + U x Jo Jn(t) 

We have proved (f!T29]h Finally we observe that from (f2725|) and (j2~26l) we obtain ([HIZO]) . □ 

Remark 2.6. In a similar way, using this time 

A+(up(«)0)(t) = A+(p(u)0)(t)«(t + r) + p(u(t))^(t)A+(u)(t) 

and A+(J(ti))(t) > p(u(t))A+(u)(t) one can prove that the opposite inequality in (|2.26p holds, and we 
have equality. Note also that equality holds also in the entropy conditions (|2.20p . 

With some additional regularity on the initial condition, one has that Uf is a Radon measure in 
(0, T) x JR. Indeed, the following proposition follows immediately from the results in |10t I21j. 

Proposition 2.7. Let uq G L°°(]R), uq(x) > k > for x G [a, b] and uq = outside [a, b]. Assume that 
uq G W 2 ' (a,b). If u is the entropy solution of (jl.3p with initial data u$, then ut is a Radon measure 
in (0,T) x M. 

From Proposition 12.71 and the results in [21], it follows that [z • z/ n W] = on d£l(t) for almost 

any i G (0, T). This permits also to define the notion of normal trace of a(v, v v ) in the sense of |2H 123] . 



3 Regularity for touching-down initial data 

Let us start by getting local estimates. 

Proposition 3.1. Let uq G L°°(M) with u${x) > n > for x G [a, b], and uq(x) = for x G" [a,b]. 
Assume that uq G W^^°(a,b). The entropy solution u(t,x) of (11.3P with u(0) = uq satisfies (i) and (ii) 
in Theorem [LTJ 

Proof. Let uo$ G L°°(iR) with uos(x) > k > for x G [a, b], uo$(x) = for x G" [a, 6], tio5 — > locally 
uniformly in (a, 6) as 5 — > + , and uo<5 G W 1,00 ([a, b]) with uniform local Lipschitz bounds in (a,b). 
Let vosiv) t> e the functions obtained by the change of variables (I2.2j) (with t = 0). Let us(t,x) be the 
entropy solution of (II. ip with U5(0) = «o5- By Theorem 1 1.1 1 we know that each us(t,x) is smooth inside 
(a, 6). Let us note that the local bounds on us and its derivatives do not depend on 5. It suffices to 
observe that this is true for the associated functions vg(t,rj) which are solutions of (I2.3p . (|2.4p . with 
initial data vg(0,ri) = vosiv)- Note that the bounds in Steps 1, 2, 3 in the proof of Theorem 12.11 are 
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independent of 5. By Remark 12.21 the Lipschitz bound in Step 5 depends only on the local Lipschitz 
bounds of vqs(v) an d are, thus, uniform in 5. Step 6 proves uniform (in 5) interior bounds for any space 
and time derivative of vs(t,rj). By passing to the limit as 8 — > 0+ we conclude that u(t,x) is smooth 
inside its support and (i) and (ii) in Theorem 11.11 hold. □ 

We now generalize our main results to initial data vanishing at the boundary of the support. 

Proposition 3.2. Let u$ G L°°(1R) with uq(x) > for x G (a, b), and uq(x) = for x (a, b). Assume 
that uq G W^°(a,b) and uq(x) — > as x — > a,b. The entropy solution u(t,x) of (|1.3|) with u(0) = uq 
satisfies (i) and (ii) in Theorem 11.11 Moreover, if uq(x) < A(b — x) a (x — a) a for some A,a > 0, then 
u(t, x) < A(t)(b + t — x) a (x — a + t) a for any x G (a — t,b + t), t > and some A(t). In that case, 
u(t, x) is a continuous function that tends to as x —> a — t,b + t. 

Proof. Let uq$ G L°°(]R) with uqs(x) = uq(x) + 5 for x G [a, 6], uo^x) = for x [a, 6], and iio<5 G 
W 1,0O ([a, 6]) with uniform local Lipschitz bounds in (a,b). Let vos(rj) be the functions obtained by the 
change of variables (|2.2|) (with t = 0). Let ^(t, x) be the entropy solution of (|1.1|) with «5(0) = uqs- By 
Theorem 1 1.1 1 we know that each us(t, x) is smooth inside (a, 6). Let us note that the local bounds on us 
and its derivatives do not depend on 5. Again, it suffices to observe that this is true for the associated 
functions v§(t,r)) which are solutions of (|2.3p . (|2.4p . with initial data v$(0,i]) = vos(r]). 

The L p bounds follow from Step 1 in the proof of Theorem 12. II for p G [1, oo) and they only depend 
on the L p bound of v $. Actually, we have 

/ v QS (rj) p drj= [ -A— T dx, 

Jo Ja uos(x)p 1 



that depends on the integrability of at the boundary points. But multiplying (|2.3p by v^cf) where 

p G [l,oo) and is a positive smooth test function with compact support in (0, 1) we obtain 



1 d 

p+ldt 



Vs +1 (t,v)4>dv+ [ \(vs) v \<l> < P [ v p s +1 4>dri+ [ v^\4> v \ drj. 
Jo Jo Jo 



Thus we derive local LP bounds for v$ which are independent of 5. We also obtain local bounds on 
the total variation of v$ which are independent of 5. To obtain a local L°° bound independent of 5 we 
observe that this follows from the identity v§(t,ri) = — h~~ ri where x = ip$(t,rj) is given by (12.20 . since 
we know that us(t,x) is locally bounded away from zero in its support [6]. Thus Steps 1, 2, 3 hold 
in their local versions. By Remark 12.21 the Lipschitz bound in Step 5 depends only on the uniform 
local bounds on v$(t,r]) and on the local Lipschitz bounds of vos(r]) and are, thus, uniform in 5. Step 
6 proves uniform (in 6) interior bounds for any space and time derivative of vs(t,rj). By passing to the 
limit as 5 — > 0+ we conclude that u(t,x) is smooth inside its support and (i) and (ii) in Theorem 11.11 
hold. The last assertion is a consequence of the comparison principle using Lemma 13.41 below. □ 

Remark 3.3. Note that the last assertion implies that if the initial profile is not vertical at the 
boundary at t = it remains non-vertical for any t > 0. Moreover, during the proof we have observed 
that if Uq has a vertical profile with ^ G L p (a, b), then u ^ x ^ G L p (a — t,b + t) for any t > 0. Thus in 
that case u(t, x) has a vertical profile at the boundary of its support. 

Due to translational invariance of (jl.3p . we state our next Lemma in an interval symmetric around 
zero. 
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Lemma 3.4. Let U(t,x) = A(t)(R(t) 2 - x 2 ) a where R{t) = R + t, a > 0. If A'(t) > 0, then U(t,x) 
is a supersolution of (|1.3I) , 

Proof. Computing the derivatives, we get 

U t = A'(R(t) 2 - x 2 ) a + 2AaR(R 2 - x 2 ) a ~\ 
U x = -2Aax(R 2 - x 2 )"- 1 , 

and 

(U 2 + U 2 x ) 1 / 2 = A{R 2 -x 2 ) a - 1 Q{x), 

1 /2 

where Q{x) = {{R 2 — x 2 ) 2 + 4a 2 x 2 ) , and then 

UU X 2Aax(R 2 - x 2 ) a 



(U 2 + U 2 ) 1 / 2 Q 

Thus, the claim 

U t > 



uu x 



Vu 2 + (u x ) 2 



holds if and only if 



jive - *r + 2Acr^ - x r- i >- 2Aa(R '-^ + ^ - ^ 



+ 



Q Q z 
AAa 2 x 2 (R 2 -x 2 )^ 1 



Q 

Let us prove that 

Indeed, the above inequality is implied by 2R > 4ax 2 /Q and 4ax 2 < (2R)2a\x\ < 2RQ. Now, we 
choose A such that 

A'(i? 2 - x 2 y > 2Aa{R2 ~ x2)a + 2Aax ( R2 ~ x2 ) a Q* 



Q Q 2 
that is, 



A , 2Aa 2AaxQ x 2Aa I xQ x \ , , 

Noticing that 

xQ x _ 4a 2 x 2 - 2x 2 {R 2 - x 2 ) Aa 2 x 2 

~Q~ Q 2 -"g^- 1, 

hence (I3.ip holds if A' > 0. We have proved that if A' > 0, then U(t,x) is a supersolution of fjl .3H . □ 
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4 Numerical experiments and heuristics 



In this section, we will propose a numerical scheme for more general equations than the RHE (11. II) . 
We deal with the Cauchy problem for the generic porous media relativistic heat equation (RHEm) [22] 
given by 

( u m u x \ 
ut = ====== (4.1 

V^ 2 + K) 2 / x 

with initial data no a probability density with compact support. In order to propose the numerical 
scheme, we make use of the change of variables to Lagrangian coordinates. As in the introduction, let 
us denote by F the distribution function associated to the probability density u and tp(t, rj) its inverse 
or generalized inverse, defined by 



<P(t,r}) :-- 



-oo rj = —\ 



inf{x: F(t,x) >r? + i}, 77 e (-|, |) (4.2) 

+00 r l = \- 



Here, we have preferred to shift the mass variable to the interval (— 2> 2) *° simplify the notations about 
boundary conditions. In this way, we simply have the relation 

F(t,ip(t,rj))= V , V e(-H). (4.3) 

For simplicity, most of the numerical tests have been chosen for even initial data. Observe that this 
change of variables is a weak diffeomorphism in case of connected compactly supported smooth u, say 
on the interval (—A(t),A(t)) in which case 

lim <p(t,rj) = ±A(t). (4.4) 
Straightforward computations show that the equation satisfied by tp in (— ^, ^) is 

1 V™" 1 ( 1 



ip t = - " ' -^ =p, (4.5) 

while at the boundary, formally, by (|4.2p and (|4.4p . we have to impose 

<A? (<. \) = (4.6) 

Moreover, thanks to the vertical contact angle property (see (|2.ip for the RHE and |22| for the RHEm), 
we have that 

lim ( — ] (t,r/)==Foo. (4.7) 



The purpose of this section is two-fold. On one hand, we heuristically observe some qualitative 
properties from the Lagrangian viewpoint. On the other hand, these properties are confirmed by 
numerical experiments with the use of an adaptation of the algorithm proposed in [16] for general 
equations in continuity form for the 2-dimensional case. 
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4.1 Numerical Method 



Equations (jl.3p and (|4.ip have been numerically treated in [341 [38] using the connection between 
nonlinear diffusions and Hamilton-Jacobi equations and numerical methods for conservation laws and 
in [10] using an appropiate WENO scheme. Here, we propose a completely different approach based on 
the optimal transportation viewpoint. As we already mentioned in the introduction an explicit Euler 
discretization of the equation satisfied by the generalized inverse (|4.5p coincides with the variational 
scheme introduced in |29[ 133]. Moreover, the theoretical result proven in [35] shows that this scheme 
applied to (|1.3p is convergent for initial data compactly supported smooth in their support and bounded 
below and above. Therefore, we plan to use a similar algorithm for Eq. (|4.ip . This Lagrangian 
formulation in ID for nonlocal and nonlinear diffusion problems was numerically analysed in [28} [T3] . 
These Lagrangian coordinates ideas were generalized to several dimensions in [16] , 

The advantages of this method are the adaptation of the mesh to the mass distribution of the 
solution in an automatic way, the immediate positivity of the solutions, and the decay of the natural 
Liapunov functional of the equations. We refer to |16] for more details and discussions on these issues. 

Here, we propose an adaptation of the algorithm in [16]. First of all, the discretization in the mass 
variable has been treated by finite difference approximations of the derivatives of the unknown ip. We 
consider a partition {r]i}i = i : N of the spatial interval [— ^, ^] and we let Aj := r^+i — r^. Note that, due 
to (|4.2p . first derivatives at the points corresponding to the nodes r/2 and i]N-i have to be taken from 
the inside of the domain. In order to avoid higher errors in the approximation of the derivative at the 
boundaries, we decide to approximate (p„ as 

( (p(m+i) - <p(th) ., , ,.s 
— if T}i < T){t) 



Vim) - <p(ra-i) 



A 



i-l 



if rji > r)(t) 



with r](t) to be specified. The derivative of the term — is computed in the other direction for better 
stability properties of the approximation of (((A?) -1 ) • At the boundary we just impose ([4. 61) . 

As explained in [16], the point r](t) has to be taken as the global maximum for u, which can be 
tracked at any time step. In all examples computed, initial data are taken to be radially symmetric and 
decreasing from the point x = rf = 0. In all of them, the global maximum stays at x = rj = 0. Therefore, 
we choose to take an even number of points N in the discretization and to take a symmetric partition 
{r]i}i = i : N of the spatial interval [—5,5]- Let us point out that the spatial partition is never uniform since 
the change to Lagrangian coordinates produces the accumulation of nodes near the global maximum. 
We instead want to follow some particular features of these type of equations such as propagation of 
fronts with a vertical contact angle or formation of singularities. Therefore, the partitions will be chosen 
accordingly in order to accumulate more points around the points ±i and other points of interest. The 
time derivative is evaluated through a simple explicit Euler scheme with the CFL condition proposed 
in [16]; i.e: 

1 



At 1 

< 



, (At?) 2 a C FL 

with acFL > 2, for the porous- medium equation which is the large-time limit behaviour of (|4.ip . see 
|22j and subsection 5.3. All our simulations are done with olcfl = 8. Although the CFL analysis 
in [16] applies only to equations written in variational form that includes (I4.ip only for m = 1, all 
numerical tests seem not to be affected by the chosen CFL condition. Finally, we point out that 
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u(t,<p(t,r]i)) = u(t,Lp(t,rjN)) = 0. Because of this fact, in all the plots which follow, the first and last 
nodes are never plotted. 



4.2 Formation of discontinuities 

4.2.1 Propagation of the support of solutions and waiting time phenomenon 

Observe that Eq. (j4.5|) and (j4.7|) imply that the speed of propagation of the support is exactly 

(\ m— 1 
/ 1T J = ±u m -\±A(t)) . (4.8) 

(here and from now on f(a ) := lim x _ s>a ± f(x) for a generic function / and point a). This coincides 
with well-known results in [21]. If we let < VK 7 ?) = ^iv) = u { L P{ r i))^ then (14.51) transforms into: 



Note that 

^(t, »7) = ( — ) (t, V) = (—) (t, 'Pit, V)) • (4-10) 

In case u(±A(t)) ^ or n x (±^4(t)) ^ if u(±^4(i)) = 0, then the boundary condition for ip is just a 
vertical contact angle using (|4.7|) - (|4.10p : 

V^,±± T ) = Too. (4.11) 




Figure 1: Left: initial datum. Right: Evolution of u$ in case m = 1 at different times. 

Consider now m = 1. By HUD, M < 1 and ift{±\ T ) = ±1, it follows that (<p v ) t (±\) > 0. This 
implies that ijjt(±7^) < by definition of ^(t, rj). In particular, this shows that in case ip(to, ±| ) = 0, 
this condition remains true for all time as shown in Proposition 13.21 

We define next w(t, rj) := ip(t, rj)il) v {t, rj) = u x (t, <p(t, rj)). The analysis above also shows that, in case 
ip(to, ±1 ) = 0, then |io(t, ±|)| < \w(0, ±|)|. On the other hand, in the bulk, w verifies the following 
equation 

1p 5 W m 3lpW 2 2 2 4 4 

W t = — 3- H 5\* w w r]W — w riW — W )■ 

(tp 2 + W 2 )2 (ip 2 + W 2 )2 
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Thus, if wq is initially bounded, w remains bounded in [—5, ^] as proved in Section [3l Observe that at 
a point 770 of maximum of w, we have 

Mvo) < -7— rrrw) ^ °) 

implying the claim. 

We show a numerical experiment with uo(x) = (1 — |cc|)+ as initial datum which does not satisfy the 
conditions of Theorem 11.11 We take N = 1000 for the simulations. We point out that since the initial 
datum is at the extremes of the support, we need a lot of nodes in the discretization near them since 
due to the change of variables (|4.3p . then ^(ii ) = =Foo and we want the numerical scheme to be able 
to capture this feature. We report in Fig. [I] the precise evolution of the support showing the smoothing 
effect at x = 0, the boundedness of the derivative all over the support including the boundaries, and 
the expansion of the boundary at precise unit speed as expected by the theory in Theorem 1 1.1 1 and the 
heuristic arguments above. 



-time=0.01 
»time=0.05 

'"'time=0.075 




Figure 2: Evolution of uq in case m = 1.5 at different times. Left: Before the discontinuity at the tip 
of the support appears. Right: Evolution of the discontinuity front after. 

Let us now take m > 1. In case u(±A(0) T ) = (i.e. ip(0,±± T ) = 0), then flO) implies that 
the support of the solution does not move at all whenever u(±A(t)) = 0. The solution will become 
positive at the tip of the support u(±A(t)) > with t > to > if and only if ipt(to, i| ) £ (0, +00] 
with u(±A(to) T ) = ip(to,±T; T ) = 0. In case u x (±A(to)) ^ 0, we can use (14.1ip to approximate terms 
(1 + V^) 1 ^ 2 — V'r? around ±5 in the expression of (I4.9P to get 

^t(*o,±i)= lim V^Oo,7?)= lim (m - l)^""^)^, ??)• 

As a consequence, ipt{to, ) > if and only if 

lim (V m+1 yto,r7)>0. (4.12) 

Observe that this condition in (|4.12p is implied by 

lim (u m+1 (x)) x (t ,±A(t)) >0. 
x-s>±A(t) 

In such case, the solution becomes positive at ±A(t) and then, according to (]4.8p . its support starts to 
increase. We note that this waiting time phenomenon is similar to that of the classical porous medium 
equation but the condition for the support to start moving is completely different to the one obtained 
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' time-0.133 
p time=0.2 



-f.5 -1 1 1.5 

Figure 3: Evolution of no in case m = 3 at different times. Top left: Before a discontinuity on the bulk 
appears. Top right: After the discontinuity front forms till it reaches the tip of the support. Bottom: 
After the discontinuity front starts to move. 

in (TTJ. Supposing a potential growth of ip, i.e. ijj(to,r]) ~ C (to, \ — |r?|) P , p > 0, for r/ — > ±^ T , then 
we obtain that ± = +oo if and only if p < 

We point out that this behavior has already been numerically obtained in [10] . In Fig. [21 we show 
this waiting time phenomenon for m = 1.5. One can observe that initially the support does not move 

since the behavior near the boundary is tp(0,7]) ~ C (0, ^ — |r?|) 5 , then the derivative at the boundary 
builds up until the behavior at the boundary reaches the critical value producing the lift-off of the 
boundary point. More interesting is the case m = 3 which we show in Fig. [3j There, a discontinuity in 
the bulk appears before the support starts to move. 



4.2.2 Formation of discontinuities in the bulk 

In view of the first example in the last section, one may think that discontinuities may appear only as 
a consequence of the waiting time phenomenon; i.e. particles tend to dissipate but their support does 
not move, which may create the discontinuities. In this section we heuristically study that it is possible 
to create discontinuities inside the bulk even if the solution are far away from zero as seen in Fig. [3l 

First we treat the case m = 1. In case of an upwards jump discontinuity or a vertical angle at a 
point 770 e] — ^,^[ such that ipnivo)^ 1 = +°o , then we also have iptivo) = — 1> Since \ip t \ < 1, then 
^Ptivo) = — 1 implies that tpt is nonincreasing to the left and nondecreasing to the right of r/o, i.e., 
((^t?) )t ^ and ((<^r?) + )t > 0. This shows that (ip(r]o)~)t > while (ip(r]o) + )t < 0, which implies that 
the size of the discontinuity reduces in for an upwards jump discontinuity or that no discontinuity is 
created if initially there is a vertical angle. 

This last phenomenon is not true if m > 1 in the case of a vertical angle at a point 770 G] — ^, 
such that iprfivo^ = +00. From the equation (|4.9|) for ip as in previous subsection, we deduce that 
ipt(Vo) = (m — l)^ m ^(r7o), and thus, a discontinuity is created. Once we have a discontinuity at r/o 
the evolution is theoretically unknown. 
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Figure 4: Evolution of solutions corresponding to uq. Top left: Initial datum uq. Top right: Evolution 
for m = 1 at small times. Bottom left: Evolution for m = 1 for larger times. Bottom right: Evolution 
for m = 4. 



In order to show this behavior we have taken two types of initial datum with N = 1000: 

3 



1 



Mx) ■■= -7X1-1,1] + 



V 2 N *[- 



- -] 

2'2> 



and 



fio(z) == i X [ _ 1 ,_|] U[ |, 1] + t 



4 X ]" 



2 ' 2 I 



We imposed a high concentration of nodes around the vertical angles or discontinuities (i.e. x 



0.8 



= initial datum 
time=0.01 
time=0.05 




-0.5 



0.5 



-0.5 



0.5 



1.5 



= time=0.01 
time=0.1 



0.8r 




-9.5 -1 -0.5 0.5 1 1.5 -1 -0.5 05 1 

Figure 5: Evolution of solutions for u$. Top: Evolution for m = 1 at different times. Bottom: Evolution 
for m = 2 at different times. 

In Fig. H]we observe the evolution of the solutions corresponding to the initial datum uq, demon- 
strating the above heuristics. In Fig. [5l we see how an initially discontinuous initial datum uq is 
smoothed during the evolution both for m = 1 (as heuristically deduced before) and for m = 2. We 
observe that the smoothing of the discontinuity is slower with m > 1 than that of m = 1 . 
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4.3 Asymptotic behavior 

In this Section, guided by heuristics, we numerically observe the asymptotic behavior of solutions to 
(|4.1j) and the rate of convergence towards their asymptotic steady state, for which no result is available 
in the literature. Performing the classical self-similar change of variables [17] that translates porous 
medium equation onto nonlinear Fokker-Planck equations given by 



v(x,t) = e*-u(e*x, k{ek - 1)) 
with k = then equation (|4.ip transforms into 

v m Vv 



vt = div xv + 



■sjv 2 + e 2t \Vv | 



(4.13) 



(4.14) 



Therefore, formally, when t — > oo solutions of (|4.14p should converge to a stationary solution of vt = 
div (xv + t) m_1 Vt)) , i.e., to a Gaussian V(x) for m = 1 or to the corresponding Barenblatt solution 
V m (x) when m > 1 given by 



V(x) 



1 _3£ 

e 2 



and 



V m {x) 



1_ 

m — 1 ^ m ~ 1 
x 



where C m is uniquely determined by the conservation of mass. In the original variables, then solutions 
should converge to the corresponding self-similar profiles obtained from V and V m via the change of 
variables (I4.13P except time translations. To be precise, the self-similar solutions are given by 



U(\x\,t) 



e « 



'Aid 



for m = 1 and U m (\x\,t) = t m +' L ( C, 



m — 1 , 2 - 



2m(m + 1) 



for m > 1, 



where is determined as above. 



■ Initial datum 
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Figure 6: Top and Bottom left: Evolution of uq in case m = 1 at different times with N = 100. Bottom 
right: log-log plot of the estimate \\u(t) — U(t)\\i with N = 1000. 
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In the following computations we have taken uq := X[_ I I]> ^ = 100- We plot the evolution of the 
initial datum for different values of m and an estimate of the difference of u — U m in the L -norm. More 
precisely, we took \\u(t) - U m (t)\\i := jjYa=i K»i>*) - U m (xi,t)\. 

Some comments are in order. First of all, in Figure [6l we note that for m = 1, while time is 
small, the numerical solution satisfies both the linear propagation of the support property, as well as 
the vertical contact angle property. However, for larger times, these two conditions are lost during the 
computation. This is due to the fact that we took a fixed number of nodes (N = 100), and as time 
increases, this number of nodes is clearly insufficient. We have observed that by increasing the number 
of nodes (for instance to N = 1000) the time in which the numerical solution is more accurate increases. 
We can also see in Figure El that, in spite of this, the numerical solution tends to a Gaussian with an 
algebraic rate of convergence that seems to be \ the one of the heat equation. However, it is exactly 
by the same reason as before that when time increases, the rate of convergence degenerates. For this 
reason, we have included in Figure [6] the .^-convergence rate with iV = 1000. 
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Figure 7: Top and Bottom left: Evolution of uq in case m = 2 at different times. Bottom right: log- log 
plot of the estimate \\u(t) — C/2 ) ] 1 1 - 

Instead, when m > 1, the support of the solution does not propagate so fast and we can observe in 
Figures [7] and [8] how the vertical contact angle property is preserved even for large times. Moreover, 
in Figures [7] and [8] we can see how the numerical solution tends to U m for m = 2 and m = 10. In 
both cases, the rate of convergence is algebraic and, numerically, it is surprisingly seen that it might 
correspond to | in the first case and to in the second one; i.e.: the same convergence rate as for the 
porous medium equation, see [17] 133], 

4.4 Convergence toward the homogeneous relativistic heat equation 

We finally show numerically how solutions to (jl.ip , converge to solutions of the homogeneous relativistic 
heat equation 

in M N x [0, T] 
Uq(x) = Uq in M N 
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10° 10 1 10 2 -O'n 5) 10»2 



Figure 8: Top: Evolution of uq in case m = 10 at different times. Bottom left: log-log plot of the 
estimate \\u(t) — Uio(t)\\\. Bottom right: zoom of the final time interval 



when the kinematic viscosity v — > +oo as already proved in [7]. In Fig. [9] we estimate the evolution in 
time of the difference in the L 1 -norm for solutions corresponding to the initial data uq = X\—l I] f° r 
different values of v with respect to the explicit solution Uh om , given by 

1 

Uhom{X,t) - - - — %[_l_ tj l +t ] 

when v — > oo. 





Figure 9: Top left: Numerical solution at t = 1 for different values of v. Top Right: Numerical solution 
at t = 100 for different values of v. Bottom: Evolution of the L 1 -difference with respect to Uhom- 
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Appendix: A primer on Entropy Solutions 



We collect in this Appendix some definitions that are needed to work with entropy solutions of flux 
limited diffusion equations. 

Note that the equation fll .3 j) can be written as 



u t = b(u, u x ) x , in Q T = (0, T)xM (A.l) 

where h(z,£) = V^/(z,^) and 

f(z,0 = zy^+W. (A.2) 

As usual, we define 

h(z,t)=b(z,t)-Z=- F M=. (A.3) 
Note that / is convex in £ and both /, h have linear growth as |£| — > oo. 



A.l Functions of bounded variation and some generalizations 

Denote by C N and % N ~ l the A-dimensional Lebesgue measure and the (N — l)-dimensional Hausdorff 
measure in ]R N , respectively. Given an open set in 1R N we denote by V(Q) the space of infinitely 
differentiable functions with compact support in £1. The space of continuous functions with compact 
support in M N will be denoted by C C (M N ). 

Recall that if f2 is an open subset of M N , a function u E L 1 (il) whose gradient Du in the sense of 
distributions is a vector valued Radon measure with finite total variation in fi is called a function of 
bounded variation. The class of such functions will be denoted by BV{Vl). For u E BV(£l), the vector 
measure Du decomposes into its absolutely continuous and singular parts Du = D ac u + D s u. Then 
D ac u = Vu C , where Vit is the Radon-Nikodym derivative of the measure Du with respect to the 
Lebesgue measure C N . We also split D s u in two parts: the jump part D^u and the Cantor part D c u. 
It is well known (see for instance pQ) that 

D j u = {u + - u-^U™- 1 L J u , 

where u + (x),u~(x) denote the upper and lower approximate limits of u at x, J u denotes the set of 
approximate jump points of u (i.e. points x E for which u + {x) > u~(x)), and v u (x) = ^jy^{x)-> 

being the Radon-Nikodym derivative of Du with respect to its total variation \Du\. For further 
information concerning functions of bounded variation we refer to [1] . 

We need to consider the following truncation functions. For a < b, let T a ^(r) := max(min(6, r), a), 
T ib = T a,b-l- We denote 

T r -={Ta,b ■ 0<a<b}, 
T+ := {Tl b : < a < b, I E M, T^ b > 0}. 

Given any function w and a,b E M we shall use the notation {w > a} = {x E M N : w(x) > a}, 
{a < w < b} = {x E M : a < w(x) < b}, and similarly for the sets {w > a}, {w < a}, {w < a}, etc. 
We need to consider the following function space 

TBV T + (M N ):={w€L 1 (M N ) + : T afi (w) - a E BV \M A ') , V T ah E T r ] ■ 

Notice that TBV^(]R N ) is closely related to the space GBV(M N ) of generalized functions of bounded 
variation introduced by E. Di Giorgi and L. Ambrosio in [JJ. Using the chain rule for BV- functions (see 
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for instance pQ), one can give a sense to Vu for a function u G TBV + (M N ) as the unique function v 
which satisfies 

VT 0|6 (ii) = vX {a<u<b y C N - a.e., V T a>b G %■ 

We refer to [1] for details. 

A. 2 Functionals defined on BV 

In order to define the notion of entropy solutions of (jA.ip and give a characterization of them, we need 
a functional calculus defined on functions whose truncations are in BV. 

Let f2 be an open subset of M N . Let g : Q x R x M N -> [0, oo[ be a Borel function such that 

C(x)|C| - D(ar) < *, C) < M'(x) + M\(\ 

for any (x, z, £) G x IRx 1R N , \z\ < R, and any i? > 0, where M is a positive constant and C, D, M' > 
are bounded Borel functions which may depend on R. Assume that C,D,M' G L l (Q). 
Following Dal Maso [25j we consider the functional: 

1Z q (u) := / g(x,u(x),Vu(x)) dx + [ g° ( x, u(x), (x) ) d\D c u\ 
Jn Jn V \Du\ J 

r ( r u +( x ) \ 
+ / / g (x, S ,v u (x))ds) dH N -\x), 

for u G 5V(f)) n L°°(Q), being u is the approximated limit of u pQ. The recession function g° of g is 
defined by 

5°(x,z,C) = lim i tg [x,z,- 

t->0+ V t 

It is convex and homogeneous of degree 1 in (. 

In case that is a bounded set, and under standard continuity and coercivity assumptions, Dal 
Maso proved in |25j that TZ g (u) is L 1 -lower semi-continuous for u G BV(Q). More recently, De Cicco, 
Fusco, and Verde [27] have obtained a very general result about the L 1 -lower semi-continuity of TZ g in 
BV{JR N ). 

Assume that g : 1R x 1R N ->■ [0, oo[ is a Borel function such that 

C\C\-D<g(z,Q<M(l + \{\) y(z,C)eM N ,\z\<R, (A.4) 

for any R > and for some constants C,D,M > which may depend on R. Observe that both 
functions /, h defined in (|A~2l . (TOl) satisfy CO)) . 
Assume that 

(<?(u(x),0) - 5 (6,0)) £L\]R N ), 

for any u G L l {M N )+. Let u G TBV+(M N ) n L°°(1R N ) and T = T a 6 — / G T+ For each G C C {M N ), 
4> > 0, we define the Radon measure g(u, DT(u)) by 



(g(u,DT{u)),<l>) := ft^CM")) + / (s(u(x),0) - s(o,0)) dx 

■/ {u<a} 



+ / 0(x) (g(«(x), 0) - 0)) dx. (A.5) 
'{«>&} 
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If (j) G C C (1R N ), we write cp = 4> + — 4> with (p + = max((j), 0), = — min(0, 0), and we define 
(g{u,DT(u)),^ := (g(u, DT(u)),(f> + ) - {g(u,DT{u))A~) ■ 

Recall that, if g(z, is continuous in (z, £), convex in £ for any z G M, and ^ G C l (M N ) + 
has compact support, then (g(u, DT(u)),(f)} is lower semi-continuous in TBV + (M N ) with respect to 
L 1 (iR Ar )-convergence |27| . This property is used to prove existence of solutions of (|A.1|) . 

We can now define the required functional calculus (see [H [5l [21] ) . 

Let us denote by V the set of Lipschitz continuous functions p : [0, +oo[— > 1R satisfying p'(s) = 
for s large enough. We write V + := {p G V : p > 0}. 

Let 5 G -P+, T G T+. We assume that u G TBV+(M N ) n L°°(1R N ) and note that 

X{ M <a}S(u) (/(«(*), 0) - /(a, 0)) , X {lt > b} S(u) (/(«(*), 0) - f(b, 0)) G L 1 (1R N ). 

Since /i(z, 0) = 0, the last assumption clearly holds also for h. We define by fs(u, DT(u)), hs(u, DT(u)) 
as the Radon measures given by (1A.5|) with /s(z,£) = S{z)f{z,C). and hs{z,Q = S(z)h(z,(): respec- 
tively. 



A. 3 The notion of of entropy solution 

Let Li(0,T,BV(M N )) be the space of weakly* measurable functions w : [0, T] — > BV(M N ) (i.e., t G 
[0, T] — > (w(t), (j)) is measurable for every <p in the predual of BV(M N )) such that J Q T ||u;(t)||s\/ dt < oo. 
Observe that, since BV(M N ) has a separable predual (see PQ), it follows easily that the map t G 
[0,T] — > \\w{t)\\BV is measurable. By Lj ocw (0,T, BV(M N )) we denote the space of weakly* measurable 
functions w : [0,T] BV(M N ) such that the map t G [0,T] -> ||w(t)|| B v is in L[ oc (]0, T[). 

Definition 4.1. Assume tfiai u G (L X (]R N ) n L°°(1R N )) + . A measurable function u :}0,T[xM N ->• 
iR is an enirop?/ solution of HT7]) in =]°> T [ xJrJV n ^ C([0, T]; L 1 (iR Ar )), r o , 6 (u(-)) — a G 
Ll cw (0,T,BV{M N )) for all0<a<b, and 

(i) u(0) = uq, and 

(ii) the following inequality is satisfied 

I I <t>h s {u,DT(u))dt+ [ [ 0h T (u,DS(u))dt 
Jo Jm n Jo Jm n 

<[ [ {j T s(u(t))fi(t)-h(u(t),Vu(t))-V(J)T(u(t))S(u(t))}dxdt, 
Jo Jr n l j 

for truncation functions S,T G T + , and any smooth function (j) of compact support, in particular 

those of the form 4>{t,x) = 4>i(t)p(x), <j)\ G V(]0,T[), p G V(1R N ), where J q (r) denotes the 

pr 

primitive of q for any function q; i.e. J q (r) := / q(s) ds 

Jo 
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